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I . INTRODUCTION 



Consider a nonrigid structure submerged in a fluid. 
Pressure in the fluid may induce response motions of the 
structure. Conversely, motion of the structure affects 
the pressure field. Such a system is called a coupled 
elasto-hydrodynarnic system. 

Many practical engineering problems belong to this 
/ 

category. A reservoir-dam system subjected to earthquake 
motion, oscillations in rocket fuel systems, a vibrating 
ship's hull, and structures subjected to blast loadings due 
to explosions are familiar elasto-hydrodynarnic problems. 

The particular elasto-hydrodynarnic system under consid- 
eration in this paper is a submarine subjected to a blast 
wave resulting from an underwater explosion. The submarine 
is modeled as a ring-reinforced, neutrally buoyant, circular 
cylinder with rigid bulkheads. It is assumed that the blast 
wave may be modeled as a plane acoustic wave with wave-front 
parallel to the longitudinal axis of the submarine. 

The interaction of an acoustic wave with an obstacle has 
been the subject of many investigations. Sette [1] first 
computed the pressure distribution on a rigid cylindrical 
surface due to a step-pulse. The response of a cylindrical 
shell to a plane step (or decaying) shock wave has been 
treated by various authors [2-5]. 

Due to the complex nature of the coupled response a 
realistic analytic formulation of this problem is intractable. 



11 



) 



However, numerical methods, particularly the finite element 
method, provide a satisfactory means of solution. 

The finite element formulation of structure-fluid 
dynamics was first introduced by Zienkiewicz, Irons an^r 
Nath [6]. They demonstrated that the coupled natural 
frequencies and mode shapes of systems could be accurately 
determined by this method. 

It is the purpose of this paper to demonstrate that the 
finite element method, as forwarded by Zienkiewicz and 
Newton [7], is applicable to the solution of the interaction 
of a blast wave with a rigid structure. Using the finite 
element formulation, a superposition theorem is developed 
which affords a method of applying the rigid cylinder results 
to a realistic, but simple, three-dimensional model of the 
submarine . 
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II. THE FINITE ELEMENT METHOD OF SOLUTION 



A. THE STRUCTURE DISCRETIZATION 

The finite element discretization of the structure 

problem is well known and is described in Zienkiewicz ' s 

\ 

text [8]. The discretization process results in an 
assembled system of equations of the form 

v 

[M] { 6 } + [C]{ 6 } + [K ] { 6 } = {R} . (1) 

In the above equation [M], [C] and [K] are respectively the 
mass, damping, and stiffness matrices of the structure 
calculated in the proper, consistent, manner, {6} is a vector 
of nodal displacements, the dots indicate time differentiation, 
and {R} is a vector of generalized nodal loads. 

It is convenient to divide the generalized nodal forces 
into two parts (after Zienkiewicz and Newton [7]) such that 

{R} = {F} + {P} . (2) 

The vector {F} is due to external forces and (P) is due to 
the fluid pressure on the interface. For the i^ h structural 
node 



P = / N ' p dS , (3) 

1 S 

where N^' is the appropriate shape function defining the 
displacement pattern in the direction normal to the boundary , 
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p is the pressure on the structure-fluid interface S, and 
the integration is over the interface. 

B. THE FLUID DISCRETIZATION 

From first principles it can be shown [6,9] that the 
wave equation 

V 2 p = 4 P > (*0 

c 

p 

where V is the Laplace operator, p is the hydrodynamic 
pressure in excess of static pressure and c the acoustic 
velocity of propagation, together with the necessary 
boundary and initial conditions are the governing relations 
for the hydrodynamic response of the fluid. 

The spatial discretization of the fluid problem is 
accomplished by a direct application of the Galerkin weighted 
residual process . If at any instant of time p is approximated 
by 



p - N.p. (summation on j = l,...,m) , (5) 

J J 

where m is the number of nodes, N. the j shape function 

J 

chosen to describe the spatial variation of p over the region 

and p . is a set of nodal pressure values which are time 
0 

dependent, then the i^* 1 weighted residual equation is 

/ N, [V 2 (N,p.) - \ (N.p,)] dR = 0 , v (6) 

R 1 J J c J J 
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where is the shape (weighting) function [8, Ch. 3] of 
the Galerkin process and R is the fluid region under 
consideration. Employing Green's theorem, the weighted 
residual equations may be transformed to 



/ 

R 



3Nj 

3x 



3N. 

— + 

3x 



3N. 3N 

_A : 

3y 9y 



3N. 3N . 

l ,i 

3 z 3z 



dR p. 



1_ 

2 



f N.N 

R 1 J 



dR p. 



3N . 

/ N, — -i 
S 1 



~ dS p_. 
3n 



= 0 



(7) 



This system of equations can be written in matrix form as 

[QHp} + [H ] {p} = {B} . (8) 

The matrices [Q] and [K] are the assembled "inertia" and 
"stiffness" matrices of the fluid. They are defined on the 
element level by: 





/ N N dR 
R -L J 



(9) 




/ 



R„ 



3N. 3N. 3N, 3N . 3N, 3N, 

± A. i. A. ±. J. “I pJ'D 

3x 3x 3y 3y 3z 3z J e 



( 10 ) 



where R g denotes the element region. 

The vector {B} results from the surface integral in (7) 
and can be written for the i^ surface node as 



b 



i 



/ 

S 



N. 



3p 

3n 



dS 



( 11 ) 
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where n is a local coordinate in the direction of the 
outward normal. By means of this vector it is possible to 
apply the boundary condition on the normal pressure gradient, 

i .e . , 



9p 

9n 





( 12 ) 



where V n is the fluid particle velocity normal to the 
boundary and p is the fluid density. It is noted that for 
a fixed boundary V n = 0 and the corresponding b^ s 0 . 

A • nonreflecting * or radiation boundary condition, 
necessary to model numerically an infinite fluid region, was 
developed by Zienkiewicz and Newton [7]. They demonstrated 
that for a plane wave normally incident upon a plane boundary 
the condition to be satisfied by p on the boundary S is 



9p 

9n 



+ 




= 0 



The corresponding boundary integral (11) is 



(13) 



c ' N i It dS ■ 



± / N lNj dS Pj 



( 1 ^) 



This leads to an additional "damping’ 1 term [D]{p} in the 
matrix equation (8), i.e.. 



[Q] (p } + [D]{p) + [H] (p } = (B) 



(15) 
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[D] is the assembled fluid "damping" matrix defined on the 
element level as 




( 16 ) 



where S g denotes the external boundary of the element. 

C. THE COUPLED STRUCTURE-FLUID PROBLEM 

The vectors {P} of the force term of equation (2) and 
(B } of equation (15) determine the coupling of the fluid 
and structure in the discretized problem. Recall from 
equation (3) that the nodal forces on the structure due to 
the pressure are 

P. = / N.' p dS = / N.'N.p. dS = [ / N.'N. dSjp. 

^ O O J J 



where [L] is the assembled matrix of element contributions 



S 



S 



or 



(P) = [L] {p} 



(17) 




(18) 



From equation (12) 



The forcing term (b) in equation (15) becomes 




V dS 



y 



or 



{B} = - [S ] (6 } 



( 20 ) 



with [S] the assembled matrix of element contributions 




Observe that 

p[L] T = [S] . (22) 

The complete elasto-hydrodynamic system has been 
formulated and can be summarized. Rewriting equation (l) 
with (2) and (17), 

[M] { 6 } + [C ] { 6 } + [K]{6 } = [L] {p } + {F } (23) 

describes the structure behavior. Similarly, equation (15) 
with (20) and (22) can be written as 



[Q] (p } + [D] {p } + [H ] (p } = - p [L ] T {6} . 



This governs the fluid response . 
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The complete system of equations can be written as 




Using general finite element techniques it is possible 
to generate and subsequently solve this unsymmetrical system 
of equations for the submarine-fluid problem under considera- 
tion. However, an alternative method based on this general 
formulation will be utilized to effect a solution using the 
following superposition theorem. 

D. SUPERPOSITION THEOREM 

It will be shown that the coupled elasto-hydrodynamic 
problem formulated above may be solved in two stages and that 
the structure response thus found is the same as in the 
complete unseparated problem. In the first stage the structure 
is constrained against responding to the loading imposed by 
the incident pressure pulse. In the second stage the structure 
and the fluid are initially at rest. External loads are 
applied to the structure to simulate the loading by the 
pressure pulse of the first stage and the structure and fluid 
responses are found. 

1 . The Complete Problem 

Equations (23) and (2*0 govern the structure and 
fluid response subject to the initial conditions 
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(«5(0)} 


= (5(0)} = 0 , 


(26) 


(P(0)} 


= (a) , (p(0) } = {b} , 


(27) 



where {a} and (b) are given vectors. It is assumed that 

{F} = 0 . ( 28 ) 

The solution of the complete problem requires that equations 
(23) and (24) with equation (28) be solved for { 6 ( t ) } and 
(p(t)}, subject to the initial conditions, equations (26) 
and (27) . 

2 . Stage On e 

The dependent variables will be (p'~^(t)} and [o^(t)J. 
Recall that in the first stage the structure is constrained 
against motion, i.e., 

{6 (1) (t)> = 0 . (29) 

For the fluid, equation (24) becomes 

[Q] (p (1) } + [D]{p (1) ) + [H]{p (1) > = 0 , (30) 

with the initial conditions 

(P (1) (0)} = (a) , (p (1) (0)} = (b) . (31) 
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Equation (23), governing the structure response, will be 
satisfied by (29) if it is required that 

{F (l) } = _ [L]{p (1) ) . (32) 

For the first stage one need only solve equation (30) 
for (p^^(t)} subject to the initial conditions of equation 
(31). {F is subsequently found from equation (32). 

3 . Stage Two 

The system is initially at rest and the structure is 
externally loaded to simulate loading by the pressure pulse. 
The dependent variables for this stage are {p '(t)} and 
U (2) (t)}. For the fluid 

CQ] (p (2) ) + [D]{p (2) ) + [H]{p (2) > = - p [L] T {6 (2) } , (33) 

with the initial conditions 

(P (2) (0)} = {p (2) (0)> = 0 . (3*0 

For the structure 

[M]{6 (2) } + [C]{6 (2) ) + [K]{6 (2) > = [L]{p (2) } + {F (2) }, (35) 

with the initial conditions 

(6 (2) (0)} = {<S (2) (0)} = 0 . ( 36 ) 
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Specify further that 



{F (2) } = - {F (1) } . (37) 

For the second stage of the solution it is necessary 
to solve equations (33) and (35) for (p^(t)} and {6^(t)> 
with equation (37) subject to the initial conditions of 
equations (3*0 and (36). 

. Superposition 

It is now asserted that 

(P(t)> = (p (1) (t)} + {p (2, (t)J , • (38) 

and 

(6(t)> = {6 (1) (t)) + (6 (2) (t)} . (39) 

In virtue of the conditions imposed on the first 
stage pressures (p^^(t)) and displacements (6^^(t)) , 
equations (29) - (32), and on the second stage pressures 
(p^ 2 ^(t)} and displacements {<5^(t)} , equations (33) - ( 37 ), 
it follows that (p(t)) and (6(t)) as given by equations (38) 
and (39) satisfy the governing equations for the complete 
problem, equations (25) - (28). This establishes the super- 
position theorem for the coupled elasto-hydrodynamic problem. 

Implementation of the superposition theorem does not 
require that the rigid body pressure distribution or the 
structure response be obtained by the finite element method. 
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It is noted that for the submarine structural model under 



consideration a three-dimensional analysis is necessary. 
However, the superposition theorem permits the rigid body 
pressure distribution to be obtained from a two-dimensional 
analysis . 
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III. WAVE PROPAGATION BY THE FINITE ELEMENT METHOD 



To evaluate numerically the interaction of a wave with 
a structure it is necessary to propagate the blast wave 
through the contiguous fluid region. Before proceeding with 
the complicated two-dimensional interaction of a blast wave 
with a rigid cylindrical structure, the characteristics of 
the one-dimensional wave propagation are studied. The 
techniques developed and utilized for the one-dimensional 
case are subsequently employed in the two-dimensional study. 

A. ONE-DIMENSIONAL WAVE PROPAGATION 

The numerical results obtained by the finite element 
method for the one-dimensional wave propagation problem are 
easily compared to well known theoretical values and thus 
provide insight into the behavioral characteristics of the 
techniques employed. 

The physical situation to be considered is a semi-infinite 
fluid strip of unit cross-sectional area. Through this region 
a pressure wave of known shape is propagating with acoustic 
velocity c. A rigid structure is located at x = 0 and the 
radiation boundary condition is applied at x = L. Equation 

(24), with { 6 ( t ) } = 0, 

[Q] (p } + [D] {p } + [H] {p } = 0 , (40) 

governs the fluid response for the case. The initial 
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conditions {p ( 0 ) } and {p ( 0 ) } are determined by the wave form 
under consideration. 

1 . Blast Wave Model 

A blast wave , which is generated by an underwater 
explosion, is characterized by an extremely rapid rise in 
pressure to a maximum value p*, followed by decreasing 
pressure. For the purpose of investigating blast wave 
propagation it would appear that a step pulse of finite 
duration would be a satisfactory wave model. However, it is 
evident that to model a moving discontinuous function by a 
finite fixed spatial discretization is not possible. Conse- 
quently, for the purpose of this investigation the blast wave 
is modeled as a ramp rise to a constant value, i.e., a ramp- 
step wave . Although the ramp-step wave is of primary interest, 
this investigation is not limited exclusively to that wave 
form and illustrative results of tests made using different 
wave forms are also given later. 

2 . Element Selection 

To study a significant portion of a blast wave, which 
is propagating at the acoustic velocity, requires a spatial 
domain of considerable extent. This is necessary to avoid 
spurious effects caused by reflections from the (artificial) 
boundaries. Accordingly, the required linear dimensions of 
the region are approximately proportional to the time interval 
during which the structure-fluid response is sought. It is 
evident that, if a fine subdivision of the one-dimensional 
fluid strip is required to represent the wave form under 
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consideration, the number of degrees of freedom in equation 
(40) may become prohibitively large when the same subdivision 
is applied to the two-dimensional region. 

Two types of finite elements, linear and cubic, are 
considered for the spatial discretization of equation (40). 

As their names imply they respectively provide for a linear 
and cubic variation of the pressure within an element. When 
propagating smooth wave forms represented by a significant 

p 

number of elements, e.g., p(x) = p* sin (irx/L) or 
p(x) = p*cos (irx/L), both the linear and cubic elements 
produce excellent numerical results. However, in choosing 
the type of finite element best suited to represent spatially 
a propagating blast wave it is necessary to consider: a) the 

general shape of the wave ; and b) two instantaneous finite 
element representations of the wave front as it propagates 
across a region. 

As established above, the general shape of the wave 
of primary interest is the ramp-step wave . For the second 
consideration, assume that the node spacing s n is equal for 
both element types. Due to the character of the linear and 
cubic elements note that three linear elements are required 
for each cubic element in the spatial representation of the 
wave front. Nov; consider a ramp wave front that at one 
instant in time is represented spatially by four nodal values, 
i.e., the ramp rise spans one cubic element. At this 
particular instant both element types are capable of an exact 
representation of the ramp wave front . For the second instant 
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in time, let the wave advance a distance equal to one node 
space. At this instant of time the linear elements continue 
to provide an exact representation of the wave front. 

However, when the wave is in this spatial position relative 
to the cubic elements, they are required to fit a cubic 
polynomial through a function which has a discontinuous first 
derivative. This results in a distortion of the wave front. 

The above qualitative discussion explains the 
observed superior performance of the linear element when 
used in the spatial discretization of equation (*I0). It also 
serves as a reminder that higher order finite elements per 
se do not guarantee superior results and may indeed lead to 
poorer results. 

3 . Time Integration 

The time integration of the governing equation (*10) 
is computed using two finite differencing techniques, the 
Newmark-B method as adapted by Chan, Cox and Benfield [10] 
and the Houbolt method [11]. Both provide a numerical method 
of determining the time dependent response of a differential 
equation of the form 

mw + dw + kw = f ( t ) , ( *1 1 ) 

where m, d and k are constants. The difference equation 
proposed by Chan, et al. [10] is 

(m + I d + Bh 2 k)w n = (2m'- (l-2B)h 2 k)w n-1 

- (m - | d + Bh 2 k)w n _ 2 + Bh 2 (f n + (i - 2)f n _ 1 + f n _ 2 ) , (*12) 
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where h Is the time step and S is a parameter which may have 
any value from 0 to 1/4. The value $ = 1/12 is chosen for 
calculations using the stability requirements established by 
Newmark [12]. 

The starting equation is 



3 

(m+ | d + 3 h 2 k)w 1 = (m + | d - (| - B)h 2 k - (J - B) w Q 

2 2 

+ (m - j-) hw 0 + Blff^ 

+ ( ( 3 ^- B) + (J- B) S3 ) h 2 f 0 , (43) 



where w^ and w^ are initial values. 

The Houbolt method [11] is developed from a consider- 
ation of a cubic curve that passes through four successive 
ordinates to obtain the difference equation 

(2m + ^ dh + kh 2 )w n = (5m + 3dh)w ^ - (4m + | dh)w 2 

+ (m + ^ dh)w n _^ + h 2 f n . (44) 

The starting procedure for the Houbolt method employed 
in the one-dimensional propagation problem is as follows : 
assume that the wave form is known at t = 0 and that the wave 
is propagating undeformed through the fluid region with the 
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acoustic velocity c. The pressure distribution p can be 
described as 

p(x,t) = f(x + ct) (45) 

for a wave moving in the negative x-direction. It is apparent 
that once the wave form is specified at t = 0 as a function 
of x the required preceding values of pressure may be calcu- 
lated directly from equation (45) and thus provide exact 
starting values for equation (44). 

It is noted that the Houbolt method is theoretically 
stable for all step sizes, while the Newmark-B method becomes 
unstable if the limiting step size established by Newmark [12] 
is exceeded. 

As an example of the results obtained by the methods 
outlined above consider a fluid strip in which the acoustic 
velocity c = 5000 ft/sec. The region of interest extends 
5000 ft. from a rigid structure located at x = 0 and is 
divided by 31 equally spaced nodes. The integration time 
step h = 0.01 sec. At time t = 0 the wave form has the shape 
and position in the region shown in Pig. 1. As the wave 
strikes the rigid structure at x = 0 it is reflected and 
propagates to the right where the radiation boundary condition 
is applied. 

After 100 time steps (t = 1.0 sec.) the wave will be 
positioned as shown in Fig. 2. The corresponding numerical 
results are shown for the Newmark-3 integration with cubic 
elements and the Houbolt integration - linear element 
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Fig. 1. Pressure-vs . -distance at t = 0 
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Pressure-vs . -distance after 100 time 



steps 
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combination. The other combinations of integration techniques 
and elements give essentially identical results. 

It is evident that to determine the optimum integra- 
tion technique it is necessary to establish criteria to judge 
the results of tests. The criteria chosen are: a) the 

pressure-time history at the structure-fluid interface 
(impulse) and b) retention of wave shape after propagation 
through the fluid. 

If the rise time of the ramp is chosen as the 
parameter necessary to describe the wave front, criteria for 
selecting the time h and the node spacing s n can be developed. 

Assume that to describe the slope of the ramp rise 
at least two linear elements are required, i.e., 

2s 

x > — -— (46) 

r — c 

Intuitively this restriction on the node spacing seems 
reasonable because it will provide that at least two spatial 
locations in the fluid strip are used to describe the slope 
of the ramp . 

It also seems reasonable to restrict the time step 
in such a way that the wave front cannot cross an element in 
one • step , i . e . , 

h < — . (47) 



31 



The above restrictions are proposed as a general 



guide in the initial selection of the nodal spacing and the 
integration time step. They are not intended to establish 
restrictions which cannot be violated. As shown below, 
satisfactory results may be obtained when the restriction on 
x^ (^6) is violated. 



effectiveness of the integration technique, numerous tests 
were made to determine which time integration technique is 
better suited for the purpose of propagating a ramp-step wave 
Consider, for example, a triangular shape wave that 
has the following characteristics: 



where x^ is the time required for the wave to decay from the 
maximum value p* to zero. The pressure-time histories at the 
rigid structure-fluid interface obtained by both integration 
methods are shown in Fig. 3. The result shown in Fig. is 
obtained by employing a cubic finite element discretization 
with the Houbolt integration technique for the same wave form 
Consider a ramp-step wave with the characteristics: 



Using the criteria established above to judge the 
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Fig. 3. Pres sure- vs .-Time at the Rigid Structure-Fluid Interface 




3*1 



Fig. Pressure-vs .-Time at the Rigid Structure-Fluid Interface 




Fig. 5. Pressure-vs . -Tice at the Rigid Structure- 
Fluid Interface (ramp-step wave). 
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The resulting pressure-time histories at the fluid-structure 
interface are shown in Fig. 5. From this and other results 
it is concluded that the Houbolt integration technique is 
better suited for the propagation of the blast wave through 
a fluid region. 

B. TWO-DIMENSIONAL WAVE PROPAGATION 

Based upon the results of the one-dimensional study, the 
linear finite element and the Houbolt integration technique 
are selected to evaluate numerically the pressure-time 
history at the rigid cylindrical structure-fluid interface 
resulting from a ramp-step wave. The geometry of the problem 
is shown in Fig. 6. 





0 



L 



Fig. 6. Two-dimensional geometry 
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The rigid cylindrical structure, of radius a, is centered 
at the origin of the x,y coordinate system. The diagram- 
matically shown ramp-step wave is propagating through the 
fluid in the negative x-direction with the acoustic velocity 
c. Due to symmetry only the upper half plane is shown and 
the rigid body boundary condition is applied along the 
x-axis. The additional necessary boundary and initial 
conditions are specified below. It is convenient to define 
a nondimensional time parameter t = ct/a, where T = 0 is the 
time of incidence of the wave at r = a, 0 = 0°. 

1 . Boundary Conditions 

The radiation boundary condition as developed earlier 
is valid only for a plane wave normally incident upon a 
plane boundary. Clearly the wave reflected from the rigid 
cylinder will not be normally incident upon the far fluid 
boundaries. It is necessary that these boundaries be 
located a sufficiently large distance from the structure to 
insure that the results at the interface are not affected. 
This investigation and Sette's [1] results both indicate that 
the transient response of the blast wave-cylindric structure 
Interaction is essentially completed in 8a/c time units after 
incidences at r = a, 0 = 0°. As a consequence the far fluid 
boundaries are located sufficiently far from the structure 
to prevent any boundary reflected waves from reaching the 
structure in the time specified. For programming convenience 
the rigid body boundary condition is applied to the boundary 
at y = L and the radiation boundary condition is applied at 
x = ±L. 
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A significant characteristic of the radiation boundary 
condition was observed during the two-dimensional study. It 
was observed that if the constant value portion of the ramp- 
step wave is applied to the boundary nodes at x = L for three 
successive time increments, then the radiation boundary 
condition continued to let the wave enter the fluid region of 
interest. In addition the interaction of the wave reflected \ 
from the structure with the boundary does not affect the 
incoming wave. This observation permitted a significant 
reduction of the fluid region under consideration, resulting 
in a proportionate reduction in the number of degrees-of- 
freedom required for the spatial discretization of equation (40). 

2 . Two-Dimensional Mesh Consideration 

In considering the two-dimensional spatial representa- 
tion (mesh) of the fluid field by the finite element method, 
it is necessary to incorporate both the physical geometry and 
the wave propagation aspects of the problem. As described in 
the one-dimensional study the integration time step is 
logically limited by the node spacing. Equation (47) governs 
this dependence. The physical dimension of the cylindric 
structure, when considered in the context of the advancing 
wave front, also restricts the nodal spacing. Clearly, if 
the node spacing s n and structure radius "a" are chosen to 
be of the same order of magnitude and the upper limit on the 
time step is utilized, the wave front could advance across 
the structure in a few time steps and the pressure-time 
history at the interface would be of suspect value. It is 
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therefore necessary to produce a mesh which affords a degree 
of fineness near the structure compatible with the integration 
time step and wave characteristics. Additionally, computer 
storage and processing time restrictions must be considered. 

After extensive experimentation the mesh shown in 
Pig. 7 was developed. Due to geometric symmetry about the 
y-axis only the region of the first quadrant is shown. This 
mesh has the following notable characteristics : 

a) Elements which form the cylindrical structure boundary 
are ’special' 5 noded, linear-parabolic elements. At 
the interface the 'special' elements represent the 
rigid boundary by a parabolic arc and provide for a 
transition to the linear elements in the fluid field. 

The interface nodes are equally spaced at angular 
intervals of 7.5° • 

b) The 'near fluid' region is characterized by a radial 
node spacing equal to the 'far fluid' node spacing and 
an angular node spacing of 15° . 

c) The transition from the cylindrical 'near fluid' region 
to the rectangular 'far fluid' region is accomplished 
with a minimum number of distorted and degenerate 
elements . 

d) The uniform rectangular 'far fluid' region is charac- 
terized by node spacing s^ = ^ a . The appropriate 
boundary conditions as discussed above are applied at 
x = ±6a and at y = 6a. 
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The overall characteristic of the mesh is element 



size uniformity. Results indicate that such uniformity in 
both spatial directions is desirable for satisfactory two- 
dimensional wave propagation. 

3 . Integration Starting Procedure 

The Houbolt integration starting procedure, equation 
(45), used for the one-dimensional case assumed a general 
solution of the wave equation (4) of the form 

p(x,t) = f (x + ct) + g(x-ct) (48) 



with g(x-ct) = 0. This solution is valid assuming that 
f ' (x) and f"(x) exist. However, the ramp-step wave has 
a discontinuous first derivative and singular higher 
derivatives. To avoid this inc onsistcncy assume that the 
ramp wave front takes the form of a cubic curve at t - 0. 



P 



+ ax' + B(x' 



(49) 



with -|2- = a + 33(x')^ 

d X 



(50) 



where x' = x - x Q . 

If p = p* and = 0 at x 1 = £ as shown in Pig. 8 

d X 

it is then possible to solve equations (49) and (50) for the 
unknown coefficients 
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Pig. 8. Cubic Wave Front 



a = 



3p* 

TT 



and 



e = -El 

P o > 

U 3 



or 



P = o’ P 1 



(1 + I 2Ll 
K 1 2 5 



1 (— ) 3 ) 
2 d ; 



( 52 ) 



By substituting into (48) and differentiating with respect 
to time 

P| = I ¥ a (1 - ( T- )2 ) ■ <53) 

*t=o 4 ^ 
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Equations (52) and (53) determine the initial conditions 
necessary for a proper formulation of the boundary value 
problem in the region of the cubic ramp wave front. 

To incorporate the initial conditions into a starting 
procedure for the Houbolt integration technique consider a 
Taylor series expansion about t = 0 

P = p 0 + P 0 t + | P 0 t2 > (5*0 



where p^ is p 



t=0 

If h is the time step 



h 2 " 

P -1 = Pq - hp o + T p 0 » 



• h 2 “ 

P +1 = P Q + hp Q + — p Q . 



s ✓ ✓ / 



The vector {p Q } is deduced from equation (iJO), i.e.. 



<P 0 > = [Q]' 1 (-[C]{p 0 ) - [H]tp 0 >) . (56) 



Thus, equations (55) with equations (56), (52) and (53) 
provide a consistent starting procedure for the two- 
dimensional propagation problem employing the Houbolt 
integration technique. 

To minimize errors in the numerical solution of the 
pressure-time history, the wave front is positioned in the 
fluid so that at t = 0 the cubic-ramp has its initial (leftmost) 
point at r = a, 9=0°. 



4. Results 



Consider the following integration, wave and mesh 
characteristics : 



h = 



_ a 



Fc 



T = 0.97a 
T r c 5 



s n = 2 a 



Implementation of boundary conditions, mesh and 
starting procedures described above results in the pressure- 
time histories at 0 = 0° and 0 = 180°, shown in Pigs. 9 and 
10. The solid line represents Sette's [1] results displaced 
along the time axis to coincide with the incidence of the 
constant portion of the cubic ramp-step wave. Additional 
curves are given for 0 = 30°, 60°, 90°, 120°, 150° in 
Appendix A. No comparison is available for the curves given 
for 120° and 150°. Tabulated values of pressure vs. time at 
each structural node are given in Appendix B. Figures 11 and 

12 provide a comparison of results when the time step is 

a o 80 a 

increased to h = and the rise time decreased to t = — 

3c r c 

The results shown in Pigs. 9 and 10 and Appendix A are 
in good agreement with those calculated by Sette [l]. This 
agreement establishes that the Houbolt [11] integration 
technique and the linear finite element are a suitable 
combination for evaluating numerically the interaction of a 
blast wave with a rigid cylinder. 




Fig. ? . Pressure- vs .-Time at the Rigid Cylinder- 
Fluid Interface (0 = 0°) 
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Fig. 10. 



Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (0 = l80°) 




Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (0 - 0°) 
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Fig. ll. 
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Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (8 = 180°) 
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Fig. 12. 



IV. THE SUBMARINE STRUCTURAL MODEL 



A simple three degree-of-freedom submarine structural 
model is developed. This model demonstrates the implementa- 
tion of the superposition theorem using the pressure-time 
history obtained from the two-dimensional finite element 
solution of the interaction of a blast wave with a rigid 
cylinder . 

A single, ring reinforced section of a submarine, 
terminated at both ends by bulkheads, is modeled. The effects 
of the stiffening rings are included by treating the shell 
as orthotropic. Each bulkhead is rigid in its own plane. 

The added mass effect due to the surrounding water is 
included by considering the water to be incompressible and 
using strip theory . 

The three modes considered are: 1) rigid body mode, 

2 ) shear beam mode and 3) buckling mode. The equations of 
motion take the form 

[M]{q} + [ K ] { q } = {f } , (57) 

where [M] and [K] are the structural mass and stiffness 
matrices respectively, {q} the vector of generalized displace- 
ments and (f) the vector of generalized forces. 

A point on the shell, Fig. 13, is located by specifying 
z and 0 . 
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Pig. 13. Submarine Structure Geometry 



The folloxving symbols are used to describe the structure 
model : 

a = shell radius, 

h = shell thickness 

s 

& = bulkhead spacing, 

u = radial displacement 

v = tangential displacement 

w = axial displacement 

E = Young's modulus of elasticity, 

G = shear modulus of elasticity 

v = Poisson's ratio. 

Values of the structural parameters used in calculations 
are specified in Appendix C. 
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The generalized modal displacements are defined as: 

1) Rigid body mode, 

u = -q^ cos 6 , v = sin 6 , w = 0 . (58) 



2) Shear beam mode, 

TT Z 

u = -q ^ sin — cos 0 , 
w = 0 . 

3) Buckling mode, 

. 2 ttz 

u = -q ^ sin -j~ cos 20 
w = A sin cos 20 



TTZ . . 

v = q^ sin — sin 0 



(59) 



v = 



1 . 2 
2 q 3 sin 



•tt sin 20 



( 60 ) 



The parameter A is evaluated by minimizing the modal strain 
energy . 

In the rigid body mode and the shear beam mode, each 
cross-section (z = constant) moves as a rigid body in the 
athwartship direction. In the buckling mode a cross-section 
is deformed, but the tangential strain = 0 . 

A. ADDED MASS EFFECTS 

Conceptually the effects of fluid pressure p m resulting 
from the radial motion of the structure should be assigned 
to the generalized force vector. It will be shown that the 
corresponding components are linear functions of q'-^ , q^ and 
q^ so that it is appropriate to include these effects by 
augmentation of the mass matrix. 
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Consider first a shell segment of infinitesmal length 



(dz) for which the radial displacement is 



u = cos e + u 2 cos 26 , 


(61) 


where u-^ = -q^ - q 2 sin , 


(62) 


, , 2 nz 

and u 2 = -q^ sin — 


(63) 


The pressure p exerted on the hull as a result 
m 


of 


radial motion is determined using potential theory, 
incompressible fluid 


For an 


<J 

ro 

xi 

B 

u 

o 

♦ 


(64) 


Let p m " p i + p 2 


(65) 


If P 1 = p 0 r COS 9 5 


(66) 


2 

and p 2 = p 0 ^2 COS 20 » 

r 


(67) 


then p^ satisfies equation (64). 

'Prom Euler's equation of motion for an inviscid 


fluid , 


9P-, 

tt — = -pu, cos 9 , 

3r r=a 1 


(68) 


p.. = p a u^ cos 6 , 

r=a 


(65) 
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where p is the fluid density. Similarly, 



r=a 



= U 2 C0S 



( 70 ) 



Therefore, the pressure is given by 



m 



r=a 



= CP 1 + P 2 3 r=a = P a ^ u i cos 9 + \ u 2 cos 20 ) 

(71) 



Let <5Wp ' be the work done by the pressure p^ during a 
virtual displacement <5u. 

Then 

2u 

<5W ' = -dz / p Su a d0 . (72) 

P o 

Substituting 

6u = 6u^ cos 6 + 6u 2 cos 20 , 



and performing the indicated integration gives 



6W p ’ 



p 1 ** *■ 

- it p a (u-^ 6u^ + ^ u 2 <$u 2 ) dz 



(72) 



If the substitutions 



6u, = -Sq-, - Sq 



•FTZ 

2 Sln X 



(73) 



and 



6u« = - 



6q^ sin 



2irz 



(7*0 
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are made and the integration over the length carried out 
we have 



<SW p = - [6q l5 6q 2 ,6q 3 ] 



m 


2m 

IT 


2m 


m 


IT 


2 


0 


0 



0 

0 

3m 

TF 




(75) 



2 

where m = piTa H (the displaced mass). The added mass matrix 
is 



[m a ] 



m 


2m 

IT 


2m 


m 


TT 


2 


0 


0 




(76) 



B. STRUCTURAL MASS MATRIX 

For a shell having mass y per unit circumferential length 
the kinetic energy in an axial length dz is taken to be 

"1 2tt t p t p 

T ' = 5 - dz / y(u + v ) a d 0 . (77) 

s d 0 

But, from equations (58) and (59), 

u = u^ cos 0 + u 2 cos 20 , ( 78 ) 

and v = -u^ sin 0 - ^ u 2 sin 20 , ( 79 ) 
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where, from equations (62) and ( 63 )} 



• • • . TTZ 

U 1 = " q l - q 2 Sln X 
u 2 = -q 3 sin 2 ~ . 



( 80 ) 

( 81 ) 



Making the appropriate substitutions and evaluating the 
resulting integrals over 0 and the length gives 



2T 



s 



m 




q l q 2 + 





( 82 ) 



where m = 2iTua£ 
s 

= shell mass (including stiffening rings). 

Adding the bulkhead mass to the rigid body mode and 
including the added-mass-matrix the total kinetic energy 
becomes 



2T 



= (m + m s + n^) ^ (m + m s ) q^ 

+ \ (m + m g ) q 2 2 + (12m + 15m g ) q^ 2 . 



( 83 ) 



For neutral buoyancy: m = m s + m^ . Let 3 = ■— , then 

equation ( 83 ) can be written as 



2T = (q} T [M] (q) , 



(8H) 
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where 



[M] = m 



'tt (1+B) 



7 (1+e) 



f d+3) 



0 

0 



0 ^j. (12+153) 



(85) 



C. STRUCTURAL STIFFNESS MATRIX 

The coefficients of the structural stiffness matrix are 
determined from energy considerations. The strain energy V 
is expressed in the form 

2V = {q} T [K] {q} . (86) 

associated with the 
are elastically uncoupled. 



(87) 



given by equations (59). 
There is no ring bending and shell bending contributions 
are negligible, so that only membrane energy need be 
considered, i.e., = 0, ~ 0 and 



Since there is no strain energy 
rigid body mode and the other modes 
the stiffness matrix takes the form 



[K] = 



0 

0 



k 



0 

22 
0 k 



0 



0 
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1. Shear Beam Mode 



The modal displacements are 
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( 88 ) 



v sr + i — t 

'z0 3z r 30 



TT . A TTZ 

j sin 0 cos -j- 



The strain energy V ^ associated with the shear beam mode is 
given by 



2V sb - { { 0 ifze 2 h s a aedz • 



(89) 



0 0 



or 



2V 



sb 



G h g a it 

2 1 



(90) 



2 . Buckling Mode 

The strain energy of the buckling mode is determined 
in two parts: the shell bending energy and the membrane 

energy V . The modal displacements are given by equations (60) 
The shell bending energy is given by 



2V, = 



l 2 tt 

f f 

0 0 



(D k 2 + 2vD k k q + D 0 k q 2 + 2 ) a d0dz 

zz ZZ0 00 Z00Z s 

(9D 



where D„ and D Q are the flexural rigidities given as 
z u 



E h 



D = 



2 (12(1 - v 2 ) 



(92) 



and 



D 0 = a D z • 



(93) 



57 



The value a is a constant of proportionality determined by 
considering the orthotropic properties of the ring reinforced 
structure. The quantity D Qz is the 'twisting' rigidity 
given by 



D 



0 z 



1 - v 



D. 



. (94) 



where effects of ring twisting are neglected. 



The curvatures < z , k 0 and k z ^ are determined from equations 



( 60 ) as: 



K = ifu 
2 8z 2 



K e 2 
a 



a 2 

[ + u ] 

90 ^ 



and 



(95) 



2 

1 9 u 



' 0 z a 9 0 9 z 



Substituting and performing the indicated integrations yields 



ttD 



2 V, = 



b 2 



[2- tt 24 3 | + 4(1 - v ) tt 2 | 



a i J 



(96) 



The membrane energy is 



£ 2 tt 



2V. 



m 



= / / (E'e ^ 2 + G Y Z 0 2 ) h g a d 0 dz , 



0 0 



(97) 
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where 



e 

z 



3w 

3z 



2ttX 
£ q 3 



cos 



2ttz 

£ 



cos 



26 



(98) 



Y 



Z0 



3v + 1 3w 
3 z r 3 6 



(- JL 
q 3 K 2£ 



2X ■ 



sin 



2ttz 

£ 



sin 



26 



(99) 



E' = E/(l - v 2 ) . 



Performing the indicated substitutions and evaluating the 
integrals gives 



ira£h 



2 V = 
m 



^ 2 X 2 



JL + 2X 
2£ a 



G] q 3 2 



( 100 ) 



The parameter X is determined by minimizing the 



membrane energy with respect to X, i.e., 
result of this minimization is 




The 



X 



-IT 



£ 

a 



Mtt‘ 



EJ_ 

G 




( 101 ) 



The required coefficients of the stiffness matrix [K] 
are determined from equations (90), (98), (100) and (101). 

D. GENERALIZED FORCES 

The generalized force vector { f } is determined by employing 
the virtual work principle. Let <5W be the work done by the 
pressure p during a virtual displacement 6q. Then 
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2tt Z „ 

6W = - / / p (|£- 6q 

n n dC lT 



9u Xri 
1 3Qo 6q 



9u 



2 3q, ~ M 3 



6q_) a dzd0 (102) 



where 



9u 



3q x ’ 3q 2 



9u and 3u 



3q. 



are determined from the modal 



displacement equations (58-60), i.e., - 



9u_ 

9q- 



= cos e, 



9u . irz n 9u 

- 55^ = Sln X cos 6 ’ - 8qJ 

the integral over the surface 



= sin^ cos 20 



Evaluating 



6VJ = [6q 1 ,6q 2 ,6q 3 ] 
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where 



f 



f 



f 



1 

2 

3 



2ir 



a.1 f p 



0 



2a& 

TT 




2tt 

/ 

0 

2 TT 

P 



cos 0 d0 , 



p cos 0 d0 , 



cos 20 d0 



(103) 



(104) 



(105) 

( 106 ) 



The generalized forces are calculated from the two- 
dimensional finite element pressure-time history. The 
resulting generalized forces are given in Fig. 14 . The 
rigid body force f^ is in close agreement with the results 
calculated by Sette [1] and Murray [ 3 ] for a step pulse. 
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Pig. 1*1. Generalized Forces-vs .-Time 
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E. GENERALIZED DISPLACEMENTS 



The Newmark-B method, equation (*12), of time integration 
was employed to evaluate the time dependent response of the 
structure. Figure 15 shows the resulting generalized 
displacements (q^E/p*a) as a function of the nondimensional 
time parameter ct/a. 

F. STRUCTURAL STRESSES 

The structural stresses are determined from the strain 
relations utilized to determine the coefficients of the 
structural stiffness matrix. 

For the shear beam mode equation (88) gives 



For the buckling mode the stress due to ring bending is 



where c^ is the distance from the centroid of the cross- 
section to the extreme fiber. The membrane tension stress 
(equation (98)) is 



s i = T z e “ GY z e * 0 r «3 sln e cos x 



(107) 





. 2 TT Z 

q, sin cos 20 



( 108 ) 



S 



3 



c 



z 



z 




2itz „ a 

cos — — cos 20 



(109) 
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Fig. 15. Generalized Displacenients-vs . -Time 
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and the membrane stress (equation (99) is 



S 4 “ T z0 Gy z0 G q 3 ^ a ^ sin % sin 26 5 ' 11 

where A is determined from equation (101). 

Table I shows the maximum nondimensional stresses S^/p* 
as a function of ct/a. 





TABLE I. 


MAXIMUM 


MODAL STRESSES 




ct 


S 1 


CM 

CO 


S 3 


s i. 


a 


P* 


P* 


p* 


p* 


0.33 


0. 0A9 


0.035 


-0.036 


0.042 


0.67 


0.720 


0.467 


-0.469 


0.553 


1.00 


3.242 


1.880 


-1.887 


2.223 


1.33 


8.561 


4.370 


- 4.386 


5.167 


1.67 


16 . 644 


7.351 


-7.377 


8.691 


2.00 


26.700 




_ i n mn 
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V. CONCLUSIONS AND RECOMMENDATIONS 



The results of this investigation confirm that the finite 
element method provides a satisfactory numerical method of 
evaluating the interaction of a plane acoustic blast wave 
and a rigid structure. Extension of the method to other two 
and three-dimensional geometries using the techniques 
employed appears to be practical. 

The submarine structural model is a simple three degree- 
of- freedom model. It was chosen to provide an illustrative 
example of the implementation of the superposition theorem 
and not as a general model of a submarine. 

It is recognized that a more sophisticated model is 
necessary to determine the response of a submarine subjected 
to blast wave loading. To improve the model it would be 
advantageous to: 

a. Include additional modes, particularly a ’breathing' 
mode to incorporate the effects of hydrostatic loading; 

b. Distribute and elastically mount the mass of the sub- 
marine not included in the shell; 

c. Remove the restrictions on the rigidity of the 
bulkheads . 

Incorporating the above improvements and nonlinear 
structural characteristic into the submarine model will 
provide for a more complete analysis of the dynamic response 
of a submarine to blast-wave loading. 
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APPENDIX A 



PRESSURE-vs .-TIME AT THE RIGID CYLINDER-FLUID 
INTERFACE FOR 0 = 30° , 60° , 90° , 120° and 150° 

The pressure-time histories at the rigid cylinder-fluid 
interface at 0 = 30°, 60°, 90°, 120° and 150° are shown in 
Figs. 16-20. The solid line represents Sette's [1] results 
displaced along the time axis to coincide with the incidence 
of the constant portion of the cubic ramp-step wave. 

Sette [1] provided no data for 0 = 120° and 0 = 150° and 
therefore no comparison is available. Data for 0=0° and 
0 = 180° is presented in Figs. 9 and 10. 

The integration, wave and mesh characteristics are: 

ch _ 1 

a. T> 3 

ct r 

1 T = °* 97 • 

fn = 1 

a 2 • 
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Fig. 16. Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (9 = 30°) 
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Fig. 17. Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (0 = 60°) 
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Fig. 18. 



Pressure-vs . -Time at the Rigid Cylinder- 
Fluid Interface (0 = 90°) 
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Pig. 19. Pressure-vs . -Tir.e at the Rigid Cylinder- 
Fluid Interface (9 = 120°) 
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Fig. 20. Pressure- vs . -Time at the Rigid Cylinder- 
Fluid Interface (0 = 150°) 
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APPENDIX B 



TABULATED PRESSURE-TIME HISTORY 
AT THE RIGID CYLINDER-FLUID INTERFACE 



The following tabulated data is the result of a two- 
dimensional finite element solution of the interaction of 
a ramp-step wave and a rigid cylindrical structure. 

The integration, wave and mesh characteristics are: 
h = a/6c. 
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97 a/c, 
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FLUID INTERFACE (continued) 
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TABLE II. PRESSURE AT THE RIGID CYLINDER- 
FLUID INTERFACE (concluded) 



APPENDIX C 



SUBMARINE STRUCTURAL PARAMETERS 



The specific structural parameters used in the calculation 
of the coefficients of equation (57) are: 



v = 0.300, 

a = hull radius = 180 in., 

H = bulkhead spacing = 480 in., 

h = shell thickness = 1.60 in., 
s ' 

8 = m /m = 0.200 

s 

The internal ring reinforcing frames are characterized 

Frame spacing = 36 . 1 in., . 

Web thickness = 1.10 in.. 

Frame depth = 9.20 in.. 

Flange width = 9-20 in.. 

Flange thickness = 1.80 in., 

a = orthotropic constant of proportionality = 127. 
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The mass matrix is 



CM] = 



9.^05 3.592 0 . 

3.592 2.821 0 . 

0. 0. 1.102 



x 10^ Ibf sec 2 /in 



The stiffness matrix is 



[K] = 



0. 

0. 

0. 



0 . 

107.9 

0 . 



0 . 

0 . 

2.236 



x 10 Ibf /in, 



77 



LIST OF REFERENCES 



1. Sette, W.J., Pressure Produced on a Non-Yielding 

Cylinder by a Step-Pulse . V/ashington D.C.: Catholic 

University of America Press, 1951. 

2. Mindlin, R.D., and H.H. Bleich, "Response of an Elastic 
Cylindrical Shell to a Transverse Step Shock Vlave." 

J. Appl. Mech. , Vol . 20, No. 2 (June, 1952), 189-195- 

3. Murray, VI. W. , "Interaction of a Spherical Acoustic 
Wave With a Beam of Circular Cross Section." Underwater 
Explosions Research Division Report 1-55, Norfolk Naval 
Shipyard, 1955- 

. Bleich, H.H. , "Dynamic Interaction Between Structures 
and Fluid." In Structural Mechanics , edited by J.N. 
Goodier and N.J. Hoff. New York: Pergamon Press, i 960 . 

Pp. 263-284. 

5. Berglund, J.W. and J.N. Klosner, "Interaction of a Ring 
Reinforced Shell and a Fluid Medium." J . Appl . Mech . , 

Vol. 35, No. 1 (March, 1968), 139-147- 

6. Zienkiewicz, O.C., B.M. Irons and B. Nath, "Natural 

Frequencies of Complex Free or Submerged Structures 
by the Finite Element Method." Proc. Conf. Vibrations 
Civil Eng., London: Butterworths , 1965, Pp - 83-93 - 

7. Zienkiewicz, O.C. and R.E. Newton, "Coupled Vibrations 
of a Structure Submerged in a Compressible Fluid." In 
Finite Element Techniques , edited by M. Sorensen. 
Stuttgart, Germany: University of Stuttgart, Pp . 359 - 380 . 

8. Zienkiewicz , O.C., The Finite Element Method in Engi- 
neering Science . New York: McGraw-Hill, 1971- 

9. Nath, Bhaskar, "Dynamics of Structure-Fluid Systems." 

In Advances in Hydroscience , Vol. 9, edited by Ven Te 
Chow, New York: Academic Press, 1973- Pp - 85-118. 

10. Chan, S.P., H.L. Cox and W.A. Benfield, "Transient 
Analysis of Forced Vibrations of Complex Structural- 
Mechanical Systems." J. Royal Aero. Soc., Vol. 66 
(July, 1966), 457-460. 

11. Houbolt, John C., "A Recurrence Matrix Solution for the 
Dynamic Response of Elastic Aircraft." J . Aero . Sci . , 

Vol. 17, No. 9, (Sept., 1950), 540-550. 



78 



12. Newmark, Nathan M. , "A Method of Computation for 

Structural Dynamics." J. Eng. Mech. , Vol. 66 (July, 

( 1959 ) 67 - 96 . 



79 



INITIAL DISTRIBUTION LIST 

No. Copies 

1. Defense Documentation Center 2 

Cameron Station 

Alexandria, Virginia 2231*1 

/ 

2. Library, Code 0212 2 

Naval Postgraduate School 

Monterey, California 939*10 

3. Department Chairman, Code 59 2 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 939*10 

*1 . Professor R. E. Newton, Code 59Ne 2 

Department of Mechanical Engineering 
Naval Postgraduate School 
Monterey, California 939*10 

5. Lt . Donald L. Atchison, USN 1 

1^708 Poway Mesa Drive 

Poway, California 9206*1 

6. Naval Ship Research and Development Center 1 

Design Applications Branch 

ATTN: Mr. L. Nash Gifford Jr. (CODE 1725) 

Bethesda, Maryland 2003*1 

7. Naval Ships System Command (SHIPS 03) 1 

Washington, D.C. 

8. Professor C.J. Garrison, Code 59Gm 1 

Department of Mechanical Engineering 

Naval Postgraduate School 
Monterey, California 939*10 



80 



i 






15342 1 



r 



Thes i s 

A85 Atchison 

" cA Finite element solu- 

t-on of the interaction 
a P'ane acoustic 
blast wave and a cy I in- 
ane structure 

•r ? < 2 ? ■ 

M ’ 

<>* JUC (U s 




Thesis 
A85 
c .1 



l r °421 



Atchison 

Finite element solu- 
tion of the interaction 
of a plane acoustic 
blast wave and a cylin- 
dric structure. 



L 




o 2768 001 91040 9 

DUDLEY KNOX LIBRARY 







